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The  AFGL  Spectral  Model 
of  the  Moist  Global  Atmosphere: 
Documentation  of  the  Baseline  Version 


I.  INTRODUCTION 

■  recent  years  there  has  been  a  widespread  shift  from  finite  difference  meth¬ 
ods  to  spectral  methods  in  the  area  of  numerical  simulation  of  the  global  atmos¬ 
phere.  This  change  has  occurred  at  both  research  centers  (for  example,  Hoskins 

1  2  3 

and  Simmons  )  and  operational  centers  (for  example,  Bourke  et  al,  Sela  ).  It  is 

due  primarily  to  the  advent  of  the  so-called  transform  method  (Eliasen  et  al,  4 

5 

Orszag  ),  which  has  permitted  spectral  methods  to  become  competitive  with  finite 
difference  methods  in  terms  of  computer  time  and  memory  requirements.  In  view 


(Received  for  publication  13  December  1982) 

1.  Hoskins,  B.J.  ,  and  Simmons,  A.  J.  (1973)  A  multi-layer  spectral  model  and 

the  semi-implicit  method,  Quart.  J.  Roy.  Meteorol.  Soc.  101:637-855. 

2.  Bourke,  W.  ,  McAvaney,  B.  ,  Puri,  K.  ,  and  Thurling,  R.  (1977)  Spectral 

methods  for  atmospheric  modeling.  Methods  in  Computational  Physics, 

Vol.  17.  B.  Adler,  Ed.,  Academic  Press,  New  Vork,  pp.  267-324. 

3.  Sela,  J.  (1980)  Spectral  modeling  at  the  National  Meteorological  Center, 

Mon.  Wea.  Rev.  108:1279-1292. 

4.  Eliasen,  E.  ,  Machenauer,  B.  ,  and  Rasmussen,  E.  (1970)  On  a  Numerical 

Method  for  Integration  of  the  Hydrodynamical  Equations  With  a  SpectraF 
Representation  of  the  Horizontal  Fields.  Inst,  of  Theor.  Meteorology, 

Univ.  of  Copenhagen,  Report  No.  2. 

5.  Orszag,  S.A,  (1970)  Transform  method  for  calculation  of  vector  coupled  sums: 

Application  to  the  spectral  form  of  the  vorticity  equation,  J.  Atmos.  Sci. 
27:890-895. 
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of  the  current  state-of-the-art  and  with  an  eye  towards  the  future  it  was  decided 
that  the  AFGL  global  model  should  be  similar  to  other  multilayer  spectral  models 
that  is,  horizontal  variations  represented  by  expansions  in  truncated  series  of 
spherical  harmonics  and  vertical  variations  represented  by  values  at  discrete 
layers.  The  coarse  resolution  model  (rhomboidal  truncation  at  wavenumber  15 
with  six  layers)  described  in  this  report  establishes  the  baseline  version  for 
future  comparisons.  It  closely  follows  the  formulation  of  the  current  operational 
model  at  the  National  Meteorological  Center  (NMC). 


2.  MODEL  EQUATIONS 
2.1  Continuous  Equations 

We  begin  with  the  equations  of  motion  for  an  ideal  gas  in  hydrostatic  balance 
surrounding  a  rotating  spherical  planet.  The  coordinate  system  is  the  familiar 
(\,<l>,o)  system  in  which  A,  <t>  are  longitude  and  latitude,  respectively,  and 


is  the  terrain  following  coordinate  suggested  by  Phillips,  B  where  p  is  the  pressure 
and  p*  is  the  surface  value  of  p.  The  horizontal  momentum  equations  are  replaced 
by  the  equations  for  absolute  vorticity  and  divergence 

*»■*  +  ?  =  t  +  fc  *  Vxy. 

D  =  V-  5C 


where  f  is  the  coriolis  parameter,  v^  is  the  horizontal  wind  vector,  and  £  is  the 
vertical  unit  vector.  These  two  equations  are  written  in  the  form  suggested  by 

7 

Bourke 


ill  =  . 
at 


a  cos  d  L 


a A  ,  ,  aB 

ax  +  cos  *  W 


+  t  •  Vx£. 


(1) 


6.  Phillips,  N.A.  (1957)  A  coordinate  system  having  some  special  advantages 

for  numerical  forecasting,  J,  Meteorol.  14:184-185. 

7.  Bourke,  W.  (1974)  A  multi-level  spectral  model.  I.  Formulation  and  hemi 

spheric  integrations,  Mon.  Wea.  Rev.  102:687-701. 
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The  continuity  equation  is 


§?■  If  ■  «••> 

where  C  =  v  •  Vq.  Vertical  integration  of  Kq.  (11)  with  the  boundary  conditions 
cr  =  0  at  a  -  0  and  a  =  1  gives  a  prognostic  for  the  surface  pressure 


and  a  diagnostic  for  b 


b  =  o(C  +  U)  -  Cu  -  Da 


(13) 


where 

a  1 

(  )°  =  J  (  )  do  and  (  )  =  J  (  )  do 

0  0 


(14) 


The  hydrostatic  and  thermodynamic  equations  are  written  in  the  forms  that 

8  1 

are  suited  to  the  Arakawa  vertical  differencing  scheme.  '  The  hydrostatic 
equation  is 


C  Tct'K 
P 


(15) 


where  is  the  specific  heat  of  dry  air  at  constant  pressure,  *  =  R/C^,  and  the 
temperature  is  given  by  T  =  T(1(a)  +  T'.  The  thermodynamic  equation  is 


3T 

at 


1 

2  A 

a  cos  q> 


9UT' 

ax 


cos  b 


avT1 

d<t> 


+  T'D 


K  . 

-  a  o 


3Tct~k 

3a 


r1(t-c-Di+f  +  hFT 

p  p 


(16) 


8.  Arakawa,  A.,  and  Mintz,  Y.  (1974)  The  UCLA  Atmospheric  Circulation 
Model,  Dept,  of  Meteorology,  University  of  California. 
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where  H  is  the  heating  rate  due  to  the  diabatic  effects  of  convection,  condensation, 
evaporation,  and  turbulent  heat  transfer,  and  parameterized  subgrid 

scale  horizontal  diffusion.  All  of  these  are  included  in  the  parameterized  physical 
processes  that  are  described  in  Section  3. 

The  forecast  model  is  completed  with  a  prognostic  equation  for  specific 
humidity,  Q, 


ay  . 
3t 


1  [dUy  ,  ,  3VQ]  ^  .  3Q  ,  0  ,  T. 

2  dT~+  cos,pW~  +  QD  ■  °  a5T  +  b  +  h*c 

jos  p  L 


(17) 


where  S  is  the  rate  of  change  of  Q  due  to  moist  convection,  condensation,  and 
evaporation,  and  ^Fq  is  the  parameterized  subgrid  scale  horizontal  diffusion. 

Once  again  these  physical  processes  are  described  in  Section  3. 

2.3  Numerical  Method* 

Because  of  the  complexity  of  the  model  equations,  realistic  solutions  can  be 
obtained  only  through  numerical  simulation.  In  this  section  we  turn  our  attention 
to  the  system  of  finite  mathematics  that  is  used  to  solve  the  equations.  For  con¬ 
venience,  we  consider  the  numerical  methods  in  three  parts:  horizontal,  vertical, 
and  time. 

2.2.1  HOHIZONTAI. 

In  the  horizontal  domain,  we  assume  that  spatial  variations  can  be  represented 
by  expanding  each  of  the  dependent  variables  in  a  truncated  series  of  surface 
spherical  harmonics,  for  example, 

M  N 

nl\,o,o,t)  r  ^  *]  y  ]  rj™  (o,  t)Y™(\,  sin  b)  ,  (18) 

m  =  -M  n  = | m | 

where  V"'(A,  sin  o)  =  P,n(sin  o)  elmA  is  the  spherical  harmonic  of  order  m  and 
n  n  r 

degree  n,  P'" (sin  o)  is  the  normalized  associated  Legendre  function  of  the  first 
kind,  and  M  and  N  are  the  truncation  wavenumbers.  The  time-  and  o-dependent 
spectral  coefficients,  q™  etc. ,  are  complex  and  reality  of  the  various  fields 
requires  q  m  =  (-1)"1  (r?^1)  where  (  )  is  the  complex  conjugate.  The  type  of 
truncation  (for  example,  rhomboidal  or  triangular)  depends  upon  the  choice  of  N. 
For  the  baseline  model,  we  will  use  a  rhomboidal  truncation,  N  =  |m|  +  M,  with 
M  =  I 
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To  obtain  the  equations  in  spectral  form,  expansions  of  the  form  Eq.  (18)  are 
substituted  into  the  continuous  equations,  and  then  the  equations  are  operated  on 
by  the  transform  operator 

it 

2 »  T 

(  )m  =  7X-  f  f  (  )Pm(sin  o)  cos  0  do  e  dx  (19) 

n  lit  J  J  n 

0  t 

“2 

for  ail  (m ,  n). 

In  the  classical  spectral  method,  it  is  necessary  to  evaluate  the  interaction 

Q 

coefficients  (for  example,  Silberman'  )  in  order  to  determine  the  spectral  co¬ 
efficients  of  the  nonlinear  terms  in  the  prognostic  equations.  These  are  integrals 
over  latitude  of  products  of  various  combinations  of  three  associated  Legendre 
functions.  Their  evaluation  on  a  computer  is  quite  time  and  memory  consuming. 

To  avoid  these  computations,  we  use  the  transform  method  that  was  developed  by 
4  5 

Eliasen  et  al  and  independently  by  Orszag.  In  practice,  the  method  consists 
of  three  distinct  steps: 

(a)  transform  the  dependent  variables  from  spectral  space  to  physical  grid 
space; 

(b)  perform  the  required  nonlinear  multiplications  in  grid  space;  and 
(e)  transform  the  nonlinear  products  from  grid  space  to  spectral  space. 

Transform  steps  (a)  and  (c)  consist  0/  a  pair  of  transforms,  a  Legendre  transform 
(spherical  harmonic  space  to  Fourier  space  and  back),  and  a  Fourier  transform 
(Fourier  space  to  grid  space  and  back).  For  the  inverse  transform  step  (a)  the 
Legendre  transform  is  accomplished  by  actually  summing  the  Legendre  series, 
while  the  Fourier  transform  is  computed  with  an  inverse  fast  Fourier  transform 
(FFT)  algorithm.  For  the  forward  transform  step  (c),  the  Fourier  transform  is 
computed  by  the  FFT  while  the  integral  of  the  Legendre  transform  is  evaluated  by 
a  Gaussian  quadrature.  The  Gaussian  quadrature  is  exact  for  our  integrands  (that 
is,  polynomials)  given  a  sufficient  number  of  points  (Gaussian  latitudes).  For 
rhomboidal  truncation,  the  transform  grid  required  to  provide  alias -free  evaluation 
of  quadratic  nonlinear  interactions  consists  of  at  least  3M  +  1  equally-spaced  long¬ 
itudes  and  at  least  (5M  +  l)/2  Gaussian  latitudes.  For  our  truncation  of  M  =  15  we 
use  48  longitudes  (7.  5°  spacing)  and  40  Gaussian  latitudes  (~4.4°  spacing). 


ft.  Silberman,  J.S.  (1954)  Planetary  waves  in  the  atmosphere,  J.  Meteorol. 
1  1:27 -34 . 


Finally,  in  presenting  the  spectral  equations,  we  will  use  the  following 
relationships 


3Y 

dA  n 


(20) 


dP1 


_  i  n  _  m  nin  /  •  .  i  ninni  —  nin 

-  cos  *  3J-  =  nen+1Pn+1  -  <n  +  l>e„  PnM  =  Hn 


(21) 


„  -n(n  +  1)Y™ 

V2Y™  =  - , - 

n  l 


(22) 


where 


2  2 

n  -  m 


4n2  ~  1 


1/2 


2 

Following  Bourke's  notation,  the  spectral  prognostic  equations  are  (excluding 
parameterized  physics) 


~  III 

on  _ 

=  -Lm(A  ,  B  ) 
3t  n  m  m 


(2  3) 


=  Lm(B  ,  -A  )  +  "fat1)-  [Em  +  *m  +  RTnqm 
3t  n  m  m  2  n  n  0Mn 

a  L 


(24) 


9Tm 


at 


B-  =  -Lm  ( (UT1)  ,(VT')  )  +  Im 
n  \  m  m  /  n 


(2  5) 


aQ 

_ n 

at 


■  -L™  (<UQ>  ,(VQ)  )  +  J"1 
n  \  m  ml  n 


(2fi) 


at 


•  cm  -  Dm 


(27) 


15 


where 


16 


We  point  out  that  because  of  the  form  of  ,  Eqs.  (32)  and  (33)  imply  that  a  direct 
spherical  harmonic  representation  in  terms  of  U^1  and  V”1  would  require  that  one 
additional  mode  be  retained  for  each  m.  Thus,  the  Legendre  series  for  and 

V  (in  terms  of  Um,  Vm)  would  be  truncated  at  N  +  1  instead  of  N. 

•!!  n  n 

finally,  we  note  that  b  is  not  considered  spectrally  since  it  is  needed  only  in 
grid  space  for  computing  the  vertical  adveetion  terms  in  the  nonlinear  interactions. 

2.2.2  VERTICAL  STRUCTURE 

In  the  vertical,  the  model  is  divided  into  a  number  of  discrete  layers.  The 
interfaces  between  layers  are  referred  to  as  levels  and  thus  the  number  of  levels 
is  one  more  than  the  number  of  layers.  All  of  the  dependent  variables  (prognostic 
and  diagnostic)  except  a  are  considered  as  layer  variables.  The  values  of  a  are 
carried  at  the  levels.  Typical  layering  and  variable  locations  are  illustrated  in 
figure  1  (note:  a  tilde  (  )  indicates  a  quantity  that  is  evaluated  at  the  levels). 

The  current  version  of  the  model  includes  six  layers  and  seven  levels  with  a  and  a 
locations  and  layer  thickness.  A,  according  to  the  values  in  Table  1,  This  struc¬ 
ture  was  established  by  defining  four  equally-spaced  a  layers  and  then  dividing  the 
lowest  and  uppermost  layers  into  two  layers  each.  This  is  done  to  provide  addi¬ 
tional  resolution  of  the  planetary  boundary  layer  and  the  tropopause  region.  Since 
the  structure  is  set  by  specifying  the  number  of  layers  and  their  spacing,  the 


LAYER  LEVEL 


Cr |  -U  1 

- t),D,T,iJ>,LI,V| 

°k-l 

*Vi  *k-i 

17,  D,  T, 

uk 
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<Tk»l  ~  ^k*l 
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aK 
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HTTTTTTTT ~ 

Figure  1.  Vertical  Layering  and  Var¬ 
iable  Location 
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Table  1.  Current  Vertical  Structure 


Level 

a 

a 

Layer 

< 

1 

0 

0.  0622 

1 

0.  15 

2 

I).  15 

0.  1985 

2 

0.  10 

3 

0.25 

0.  3699 

3 

0.25 

0.  50 

0.  6220 

4 

0.25 

0.75 

0.  8242 

5 

0.  15 

0.  90 

0.  9497 

6 

0.  10 

| 

1.  00 

values  of  a  at  the  levels  are  immediately  prescribed.  Determining  the  values  of 
cr  in  the  layers  (that  is.  the  point  at  which  the  dependent  variables  should  be  loca¬ 
ted)  is  not  so  straightforward.  Many  researchers  have  chosen  to  simply  pick  the 
midpoint  of  each  layer  (for  example,  Bourke  ).  However,  since  we  have  chosen  to 

3 

construct  our  baseline  model  according  to  the  NMC  formulation,  we  follow  Sela 
and  determine  ct  in  the  layers  from 


CTk 


(1  +  K)(5k+1  '  ^ 


1  da1^ 
1  +  k  do 


(34) 


In  order  to  evaluate  vertical  advection  terms,  they  are  first  rewritten  in  flux 
form,  for  example. 


.  au  _  9oU  tI  ad 

°  a o  3ct  3ct 


(35) 
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The  finite  difference  analog  of  the  flux  form  for  layer  k  is 


3bU  ii  8b 

T5~  '  u  JS 


U, 


k+ 1  "k+1  k  k 


-  u. 


-  U, 


k+1 


-  bk 


(36) 


Since  U  is  a  layer  variable,  we  must  interpolate  the  layer  values  to  determine  U 

Q 

(that  is,  the  level  values).  Following  Arakawa  we  use  the  relationship 


Vi  =  V 


(37) 


and  upon  substituting  into  Eq.  (36)  we  have 


8bU  ,,  db_  _  1  r  r- 

8a  8a  ~  TS^  [ak+l 


‘V+i'V  +  W-'W 


(38) 


which  is  a  quadratic  conserving  differencing  scheme.  Similar  expressions  can  be 
written  for  vertical  advection  terms  involving  V  and  Q. 

For  the  vertical  advection  term  in  the  thermodynamic  equation,  Eq.  (29),  and 
for  the  hydrostatic  equation,  Eq.  (31),  we  apply  Eq.  (37)  to  the  entire  quantity 
Ta  *  so  that  "f  is  determined  from 


Thus,  the  vertical  advection  term  in  the  thermodynamic  equation  (in  flux  form)  is 
given  by 


k  8aTa  K 

-  8b 

_  1 

A.  j 

K  \ 

(  °l  \1 

°  8a 

1  8o 

°k+l  j 

k  k+1  k  I  +  a. 

^°k+l  /  k  ' 

\k  °k+l  Tk_1j 

(40) 

The  finite  difference  analogs  of  the  vertical  integrals  of  Eq.  (14)  are 


k+1 


n<Vl  * / 

o  j=i 


(41) 
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and 


o./<  >*,  =  £<  ,.d.  . 


We  conclude  this  subsection  by  considering  the  solution  of  the  hydrostatic 
equation.  With  the  aid  of  Eq.  (39),  the  finite  difference  form  of  the  hydrostatic 
equation  is 


4>,  -  =  o,.,  T,  .  ,  +  a,  T, 

k  k+1  k+1  k+1  k 


where 


ak+i  =  T-  1 '  y — 


ft.  -St  fe-l 


^  t  n? 


The  geopotential  of  the  lowest  layer  is  determined  from  the  integral  constraint 


J  (RT  -  (*  -  «*)]  do  =  0  , 


where  **  is  the  surface  geopotential  (that  is,  terrain  height).  In  finite  difference 
form  Eq.  (44)  is 


£  ,RTk  -  <*k  •  ■  °  • 


Upon  combining  Eqs.  (43)  and  (45)  we  have  an  expression  for  the  height  of  the 
lowest  layer, 
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htk^k  + 


(4(1) 


”K  '  *' 


K-l 

E 

k-  ! 


K-l 

nTk-Z‘v1V.t'W 

J=k 


l’lius,  the  geopotential  for  all  layers  can  be  determined  from  the  hydrostatic  equa¬ 
tion  by  solving  l.q.  (4i>)  for  the  lowest  layer  followed  by  Eq.  (43)  for  all  other 

layers . 

3.  3  mu:  integration 

l  lie  time  integration  is  accomplished  by  the  so-called  semi-implicit  method, 
which  is  a  mixture  of  implicit  and  explicit  techniques.  In  this  method,  linear 
terms  associated  with  gravity  waves  in  the  time-dependent  equations  are  treated 
implicitly  while  other  terms  are  treated  explicitly.  Such  an  approach  avoids  the 
problem  of  having  to  deal  with  a  nonlinear  algebraic  system  associated  with 
implicit  methods  for  nonlinear  differential  equations;  and  yet  it  gives  us  a  time 
integration  scheme  in  which  relatively  large  time  increments  can  be  used 
(Hubert10).  Had  we  used  an  explicit  scheme,  computational  stability  considera¬ 
tions  would  have  dictated  the  use  of  a  much  smaller  time  increment  (for  example, 

1°  min  for  a  leap  frog  scheme  as  compared  to  60  min  for  the  semi-implicit  scheme) 

There  are  various  ways  to  implement  the  semi-implicit  method.  The  one 
adopted  here  is  due  to  Bourkc.  ‘  The  key  point  is  to  eliminate  T  and  q  among  the 
discrete  forms  of  Kqs.  (24),  (25),  and  (27),  leaving  1>  as  the  sole  unknown.  The 
mathematical  details  of  the  derivation  of  the  algebraic  divergence  equation  are 
given  in  Appendix  B.  Here  we  only  outline  the  procedure; 

(a)  Express  T  in  Kq.  (24)  in  terms  of  T,  using  the  hydrostatic  equations, 

Kqs.  (43)  and  (45)  and 

(b)  Approximate  time  derivatives  in  Kqs.  (23)  through  (27)  by 

<*F  F+  -  F* 

3t  TSl  ’  ' 


Here  F  is  any  of  the  prognostic  variables,  and  F+  and  F  are  values  of  the  vari¬ 
ables  at  the  next  and  previous  time-steps,  respectively. 

(c)  Kcplace  F  in  the  linear  terms  on  the  right-hand  side  of  Eqs.  (24),  (2  5), 
and  (27 )  by 


10.  Kobert,  A.  (197!))  The  semi-implicit  method,  Numerical  Methods  Used  in 

Atmospheric  Models,  Vol.  2.  GARP  Publication  Series,  No.  fE  Geneva, 
WMO,  pp,  419-43!). 


'.S-tjet-S-  x  . . 
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F  =  <>F+  +  (1  -  ••  >F ■'  , 

where  1/2  So  s  l, 

(d)  Eliminate  T  ,  q  in  the  discrete  divergence  equation  by  using  the  discrete 
forms  of  Fqs.  (2  5)  and  (27),  thus  leaving  an  algebraic  system  in  L>+.  Once  D+ 
values  are  found,  Eqs.  (2  5)  and  (27)  can  immediately  be  solved  for  T+  and  q+. 

Since  there  are  no  linear  terms  on  the  right-hand  side  of  the  vortieity  and  moisture 
equations,  these  two  equations  are  integrated  explicitly. 

In  the  current  version  of  our  model,  the  following  time-stepping  sequence  has 
been  implemented: 

(a)  Advance  n  and  O  using  an  explicit  leap  frog  scheme, 

(b)  Advance  D  to  the  new  time  step  t  +  At  by  solving  the  linear  algebraic 
system  Eq.  (Bll)  from  Appendix  B, 

(c)  Advance  T  and  q  using  Eqs.  (B9)  and  (BlO)  from  Appendix  B, 

(d)  Apply  the  moisture  physics  and  adjust  the  predicted  values  of  T  and  Q, 

and 

(e)  Time-smooth  all  of  the  prognostic  variables  at  current  time  t  to  prevent 
decoupling  between  the  even  and  odd  time-step  values.  Thus,  the  latest  predicted 
value  of  a  variable,  F  ,  is  used  to  compute  a  time-smoothed  value  for  the  next 
time -stepping  cycle: 

F "  =  F  +  )>(F”  -  2 F  +  F+)  , 

where  F  becomes  the  time-smoothed  F  value  for  the  next  time-cycle,  and  fi  is 
a  given  constant.  For  the  very  first  time -step,  a  forward  scheme  is  used. 

For  the  results  reported  below  we  use  the  values  At  =  60  min,  a  -  1  (back¬ 
ward  implicit),  and  (J  =  0.  02. 

3.  PARAMETERIZED  PHYSICAL  PROCESSES 

To  establish  our  baseline  model,  we  have  adapted  the  various  parameterized 

3 

physical  processes  that  appear  in  the  current  NMC  spectral  model.  These  fall 
into  three  broad  categories:  boundary  layer  processes;  moisture  physics,  includ¬ 
ing  convective  adjustment  as  well  as  the  so-called  large-scale  saturation  process; 
and  subgrid  scale  diffusion. 

3.1  Boundary  Layer  Physics 

Boundary  layer  processes  are  included  in  the  model  to  simulate  the  interaction 
between  the  surface  of  the  earth  and  the  atmosphere.  In  particular,  three 
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mechanisms  are  considered  in  terms  of  flaxes  into  the  atmosphere:  surface  fric¬ 
tion  (that  is,  momentum  flux),  sensible  heat  flux  from  the  surface,  and  evapora¬ 
tion  from  the  surface  (that  is,  moisture  flux). 

7 

The  surface  friction  is  formulated  following  Bourke  with  some  minor  modifi¬ 
cations.  It  is  simulated  by  including  the  stress  terms 


JL 

P* 


% 


and 


ZJL 

p* 


(47) 


in  Eqs.  (3)  and  (4),  respectively.  For  purposes  of  the  vertical  finite  difference 
scheme,  and  are  defined  at  the  levels  and  are  set  to  zero  except  at  the 
surface  where  they  are  given  by  the  bulk  aerodynamic  formula 


.tv  =  pcd  yii  V 


where  p  is  the  density.  C  .  is  the  drag  coefficient,  E  is  the  kinetic  energy,  and  V* 
is  the  horizontal  vector  whose  components  are  the  pseudo-velocities  U,  V  at  the 
surface  of  the  earth.  Strictly  speaking,  p  and  E  should  also  be  evaluated  at  the 
surface.  However,  in  view  of  the  crude  nature  of  the  stress  formulation  and  for 
simplicity,  p  and  E  are  assigned  the  values  from  the  lowest  model  layer.  is 
given  by  V  at  the  lowest  layer  rotated  by  an  amount  a  towards  low  pressure.  In 
addition,  a  latitudinal  weighting  factor  of  sin  <p  is  included  with  sin  a  to  reduce  the 
influence  of  turning  near  the  equator  and  to  provide  the  proper  sense  of  rotation  in 
the  Southern  Hemisphere.  Thus,  the  surface  stress  terms  are  given  by 


(r,Jx  -  »K  Cd  <UK  cos  u  Vj.  sin  a  sin  </>) 

=  PK  Cd  ^  cos  o  +  Uj^  sin  a  sin  b)  ,  (49) 


where  a  subscript  K  indicates  the  value  at  the  lowest  model  layer.  The  value  of 
o  =  20  is  fixed  for  all  locations.  Values  of  C'd  vary  geographically  from  a  min¬ 
imum  of  0.  0013  over  the  open  oceans  to  a  maximum  of  0.  OOP  over  the  Himalayas. 

Sensible  heat  transfer  from  the  surface  is  included  as  part  of  the  diabatic 
heating  term,  ~~  ,  in  the  thermodynamic  Eq.  (IG>.  The  parameterization  allows 
P 

only  an  upward  heat  flux  into  the  lowest  model  layer  and  only  over  the  oceans.  The 
exact  form  of  the  heating  term  is 


H 

r~  = 

p 


JL 

P:- 


(50) 
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It 


where  F  is  defined  at  the  levels  and  is  set  to  zero  except  at  the  surface  where  it 
0 

is  given  by  the  bulk  aerodynamic  formula 

(F^)#  =  Pp,  C^,  \/2E ( Q ^  ,  (>l) 

where  Q  is  the  potential  temperature,  a  subscript  K  indicates  the  value  at  the  low  - 
est  model  layer,  and  C^,  is  the  thermal  drag  coefficient  given  by 

CT  =  cd  +  7  x  i°'5  .  (52) 

This  formulation  of  allows  the  heat  transfer  process  to  be  highly  wind-speed 
selective  (that  is,  it  becomes  more  efficient  as  the  wind  speed  increases).  For 
these  computations  we  use  p*  as  the  reference  pressure  for  defining  potential 
temperature  and  thus  Eq.  (51)  becomes 

(F^)J}.  =  Pj,  C^,  V'2Ej-.  (T#  -  T^  )  ,  (53) 

where  T*  is  the  sea  surface  temperature. 

Evaporation  from  the  surface  is  included  in  the  moisture  source/sink  term,  S, 
in  Eq.  (17).  Evaporation  is  allowed  only  over  the  oceans  and  only  if  the  lowest 
model  layer  is  not  saturated  relative  to  the  ocean  surface  and/or  there  is  no  pre¬ 
cipitation  in  the  lowest  layer.  The  mathematical  formulation  of  the  evaporation 
physics  closely  resembles  the  sensible  heat  transfer.  The  form  of  the  evaporation 
term  is 


S  =  JL 


x  Q 

p,;,  15 


where  Fq  is  defined  at  the  levels  and  is  set  to  zero  except  at  the  surface  where  it 
is  given  by 


<V*"PKCQ'/5^«s<V‘QK>  * 


where  Qs(T.,.)  is  the  saturation  specific  humidity  at  the  sea  surface  temperature, 
and  the  moisture  drag  coefficient  is  given  by 


CQ  =  CT  =  Crf  +  7  X  10  J  v  2Ek  . 


Finally,  to  provide  additional  computational  stability,  the  surface  drag,  the  heat 
flux,  and  the  moisture  flux  are  lagged  in  time. 

3.2  Moisture  Physics 

In  this  section  we  describe  the  parameterization  of  the  processes  that  are 
associated  with  the  phase  changes  of  water.  These  are  adapted  directly  from  the 
NMC  spectral  model'’  with  only  minor  modifications  to  account  for  the  difference 
between  the  sigma  structures  of  the  two  models.  The  values  of  pressure-dependent 
parameters  that  appear  in  the  algorithms  are  defined  by  simple  interpolations  of 
those  employed  by  NMC.  The  procedure  consists  of  a  sequence  of  three  adjust¬ 
ments  of  temperature  and  specific  humidity.  Each  step  represents  a  particular 
aspect  of  physics  involving  the  vertical  stratification  of  the  variables.  The  ad¬ 
justments  are  applied  after  the  provisional  values  of  the  dynamic  variables  at  a 
new  time  step  are  obtained  through  time  integration. 

The  current  practice  makes  no  distinction  between  liquid  and  solid  phases  and 
all  values  of  the  relevant  physical  parameters  such  as  specific  heat  and  latent 
heat  refer  only  to  either  gaseous  or  liquid  state. 

3.2.  1  CONVECTIVE  ADJUSTMENT 

The  first  adjustment  is  based  on  Kuo's11  simulation  of  organiz  'd  cumulus 
convection.  This  moist  convective  adjustment  scheme  consists  of  two  parts.  The 
first  part  of  the  scheme,  which  has  been  described  in  great  detail  by  Phillips,  1_ 
accounts  for  the  changes  that  occur  in  temperature  and  specific  humidity  when 
convection  redistributes  heat  and  moisture  in  a  moist  unstable  column.  In  imple¬ 
menting  the  scheme,  the  only  change  we  make  in  Phillips'  description  is  the  criti¬ 
cal  value  of  the  moisture  convergence  above  which  the  adjustment  is  invoked.  We 
-7  - 1 

use  a  value  of  1  x  10  •  At  cb  s  where  At  is  the  time  step  in  seconds.  This 

is  the  appropriate  value  for  our  choice  of  vertical  structure. 

The  second  part  of  the  convective  adjustment  scheme  takes  into  account  the 
effects  of  evaporation  of  falling  water.  This  will  cause  c  hanges  in  temperature, 
specific  humidity,  and  the  amount  of  rainfall.  Phillips  only  gives  a  brief  descrip¬ 
tion  of  this  part  and  thus,  for  the  sake  of  completeness,  wc  provide  the  balance. 

First,  we  define  the  moisture  deficit  for  each  unsaturatod  layer  in  the  unstable 
column 

11.  Kuo,  if.  I ..  (infifi)  On  formation  and  intensification  of  tropical  cyclones 
through  latent  heat  release  by  cumulus  convection,  J.  Atmos.  Sci. 
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Phillips,  N.A.  (U'79)  The  Nested  Grid  Model,  NOAA  Technical  Report 
NWS -2 2,  7  0pp. 
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(57) 


w ,  ,  =  (Q  .  .  -  y,  .  )au,  , 

uel ,  k  Stit,  k  0,  k  k 

where  Wsat  j.  is  the  effective  saturation  specific  humidity  (chosen  to  be  eight-tenths 
of  the  theoretical  value),  Q()  ^  is  the  latest  provisional  value  of  specific  humidity, 
anil  is  the  thickness  of  the  layer. 

The  amount  of  total  condensate,  K,  starts  in  the  highest  saturated  layer  of  the 
unstable  column  with  the  value  given  by 


1  k  =  1)TKIT(\  ^'k  • 


in  which  DTK UO  denotes  the  temperature  adjustment  to  that  layer  brought  about 
by  the  Kuo  convection.  As  the  condensate  falls  it  will  either  undergo  evaporation 
in  unsaturated  layers  or  accumulate  additional  condensate  in  saturated  layers. 

In  an  unsaturated  layer,  evaporation  must  continue  until  all  of  the  condensate 
evaporates  or  until  the  layer  becomes  saturated.  In  the  former  case  k>  1-V, 

and  the  adjusted  values  of  specific  humidity,  temperature,  and  condensate  are 


^Kk  =  %l,k  +  “Vl^k 


c  If.  , 
T  =  T  -  il  ^-l 
l.k  li.k  E  3o7~ 
v  k 


Fk  ’  ° 


In  the  later  case  ^  k  l’k.  and  only  a  part  of  condensate  evaporates,  thus 
bringing  the  layer  to  the  level  of  effective  saturation,  with  the  results 


^l,k  ^sat,  k 


ri,k  =  r0,k  *  Qdef,  k 

V 


K  =  K  -  UCi 
k  k- 1  Ao, 
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In  a  saturated  layer,  Q  and  T  remain  unchanged  and  F  accumulates  the  addi¬ 
tional  condensate  so  that 


(fil) 


<R2) 


where  p  is  the  surface  pressure. 

3.2.2  LARGE-SCALE  SATURATION 

The  second  adjustment,  which  is  also  described  in  Phillips'  report,  is  re¬ 
ferred  to  as  the  large-scale  saturation  process.  It  is  purported  to  adjust  tem¬ 
perature  and  specific  humidity  according  to  the  wet-bulb  relation.  There  is,  how¬ 
ever,  built-in  empiricism  in  the  formulas  currently  used  to  calculate  the  amounts 
of  assumed  and  allowable  evaporations.  A  brief  account  of  these  aspects  is  given 
here. 

For  each  layer  we  define 


where  and  Qpff  will  be  referred  to  as  the  apparent  saturation  specific  humid¬ 

ity  and  the  effective  saturation  specific  humidity,  respectively.  C  ^  ancl  ^eff  are' 

empirical  constants  that  are  determined  as  follows.  C  takes  the  value  of  0.  95 

app 

in  the  topmost  of  the  moisture-carrying  layers  and  0.  80  elsewhere.  On  the  other 
hand,  C'e^  assumes  the  value  of  unity  in  the  bottom  layer  and  the  value  0.  85  in  the 
fourth  lowest  layer  and  above.  In  the  second  and  third  lowest  layers,  however,  it 
is  assumed  that  C  jj  increases  with  temperature  quadratically  to  the  extent  that 
its  value  lies  in  the  range  between  0.85  and  0.95.  The  relationship  between 
and  temperature  is  given  by 


Ceff  =  0-85  +  x(1  +  x)/200 


ceff  S  °-95 


(64) 


1  ,  =  t  Q  . 

app,  k  app  s,k 

j 

i 

eff .  k  eff  ws,  k  ' 

(63)  * 
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where 


..  _  T  -  261.  16 

X  O 


Next,  in  each  of  the  moisture -carrying  layers,  a  deficit,  S^,  is  defined  by 


Sk  g  (Q0,k  '  Qeff,k)  Aok 


(65) 


where  Q0  ,  is  the  latest  provisional  value  of  specific  humidity.  It  is  assumed  that 
if  <  0  evaporation  will  occur  and  no  additional  precipitation  is  contributed  by 
the  layer.  The  amount  of  evaporation  is  given  by 


E 


(Q 


0,  k 


Qapp,  k*  Aak 


(66) 


in  which  C  is  an  efficiency  factor  (0  <  C  <  1)  that  decreases  downward.  The 
amount  of  precipitation  reaching  the  adjacent  layer  below  is  then  reduced  by  the 
amount  E.  If  g.  0,  on  the  other  hand,  there  is  no  evaporation  but  an  additional 
amount  of  precipitation  represented  by  is  contributed  by  the  layer 

Finally,  the  adjustments  on  temperature  and  specific  humidity  are  given 
respectively  by 


Ql,k  =  Q0.k 


T 


1,  k 


(67) 


3.2.3  ADIABATIC  ADJUSTMENT 

The  third  adjustment  insures  that  the  thermal  stratification  is  statically 
stable.  When  instability  results  from  the  last  provisional  values  of  temperature, 
it  is  assumed  that  complete  mixing  takes  place  in  both  heat  and  moisture  so  that 
the  adjusted  potential  temperature  and  specific  humidity  assume  uniform  values 
that  are  the  mass-weighted  averages  of  the  respective  quantities  within  the  un¬ 
stable  layer. 

In  addition  to  these  adjustments  that  are  based  on  physical  reasoning,  there 
are  two  other  adjustments  made  on  specific  humidity  that  owe  their  existence  to 
the  computational  aspects  of  the  model.  When  the  initial  grid-point  values  of 
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specific  humidity  are  transformed  into  spectral  amplitudes,  the  latter  are  uni¬ 
formly  reduced  to  nine-tenths  of  the  original  values.  The  other  adjustment  is 
introduced  when  the  spectral  coefficients  are  transformed  into  grid-point  values 
at  each  time  step.  Here,  whenever  a  grid-point  value  turns  out  to  be  negative  it 
is  simply  set  to  zero. 

3.3  Siiligrki  Stale  Diffusion 

The  parameterized  subgrid  scale  diffusion  in  the  model  serves  two  purposes, 
both  related  to  dissipative  effects.  From  the  point  of  view  of  physics,  this  proc¬ 
ess  is  included  to  simulate  the  assumed  dissipative  effect  of  the  unresolved  scales 
of  motion.  For  computational  purposes,  it  is  included  to  prevent  the  spurious 
accumulation  of  energy  in  the  smallest  resolved  scales  of  a  numerical  model  with 
a  finite  spatial  resolution. 

In  the  current  model,  the  diffusion  is  parameterized  by  adding  a  term  of  the 

form  -rv^(  )  to  the  prognostic  equations  for  rj,  D,  T,  and  Q.  The  v4  operator  is 

more  scale  selective  and  more  effective  in  controlling  the  higher  wave  numbers 

than  the  traditional  v“  formulation.  In  practice,  two  different  diffusion  coefficients 
17  4  - 1 

are  used:  =  4  X  10  m  sec  for  diffusion  applied  to  all  modes  of  the  diver¬ 

gence,  and  i  y  =  1  X  X01(i  m4  sec  1  for  diffusion  applied  only  to  the  upper  half  of 
the  rhomboid  (n  >  M)  for  q,  T,  and  Q.  For  computational  stability,  the  diffusion 
is  lagged  in  time.  There  is  no  vertical  diffusion  included  in  the  present  formula¬ 
tion. 


4.  DATA  PKOCKSSIM, 

Before  the  model  forecast  can  begin,  the  initial  data  must  be  prepared  in  a 
form  appropriate  for  input  into  the  model.  Similarly,  at  the  end  of  a  forecast,  the 
forecasted  data  must  be  put  into  a  form  appropriate  for  certain  standard  graphical 
and  numerical  displays. 

In  general,  preparation  of  the  initial  data  consists  of  three  steps:  (a)  objective 
analysis  in  which  the  observed  meteorological  data  are  analyzed  and  represented 
at  the  mandatory  pressure  levels  on  a  regular  latitude-longitude  grid;  (b)  pre¬ 
processing  in  which  the  analyzed  data  are  transformed  from  the  analysis  grid  at 
the  mandatory  pressure  levels  to  the  model  grid  in  the  model  o  layers  and  then 
transformed  into  spectral  coefficients  in  the  a  layers;  and  (c)  initialization  in 
which  the  preprocessed  data  are  put  through  a  filtering  procedure  in  which  certain 
undesirable  gravity  modes  (that  is,  meteorological  noise)  and  their  tendencies  are 
removed.  For  the  current  baseline  model  we  use  the  FGGK  level  111 -A.  data  that 
have  already  been  analyzed;  and  thus  we  skip  step  (a).  The  preprocessing  and  the 


normal-mode  initialization  have  been  adapted  from  the  NMC  spectral  model.  3  The 
only  changes  made  in  these  algorithms  are  those  necessary  to  reflect  the  differ¬ 
ences  in  resolution  and  a  structure  between  the  two  models. 

The  postprocessing  of  the  forecast  data  is  essentially  the  inverse  of  the  pre¬ 
processing,  that  is,  the  forecasted  spectral  coefficients  in  the  model  a  layers  are 

transformed  to  physical  values  at  the  mandatory  pressure  levels  on  a  regular 

3 

latitude -longitude  grid.  This  procedure  is  also  adapted  from  NMC. 

4.1  Preprocessing/Postprocessing 

Since  the  preprocessing  and  postprocessing  are  essentially  the  inverse  of  one 
another,  there  are  certain  features  that  are  common  to  both.  The  most  prominent 
among  them  are  the  following: 

(a)  Between  two  levels  and  at  which  data  are  available,  temperature, 
wind,  and  specific  humidity  are  assumed  to  be  piece-wise  linear  functions  of 
Z  =  In  p,  that  is, 

F(Z)  =  F(Z  )  +  A(Z  -  Z  )  .  Z  <  Z  <  Zu  ,  (G8) 

a  a  a  —  =  b 


where 


A  = 


F(Z,  )  -  F(Z  ) 
b  a 

Zu  -  Z 
b  a 


(b)  Geopotential  and  temperature  are  related  to  each  other  by  the  hydrostatic 
equation 


It  then  follows  from  (a)  that  $  is  a  piece-wise  quadratic  function  of  In  p  and  may 
be  expressed  as 


4>(Z)  =  4>(Z)  +  A(Z  -  Z)  +  5.  (Z  -  Z)2  ,  Za  <  Z  <  Zfa  ,  (70) 

where  the  constants  A  and  B  are  determined  from  the  observed  values  4>(Za), 
$(2^),  T(Za>,  and  TfZ^)  by  the  following  relations: 
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(71) 


4.1.1  PREPROCESSING 

The  first  step  of  the  preprocessing  consists  of  transforming  the  analyzed  data 
from  the  analysis  grid  to  the  model's  computational  grid.  The  FGGE  level  II1-A 
data  are  available  on  a  regular  latitude -longitude  grid  with  a  spacing  of  2.  5°  (that 
is,  144  longitudes  by  73  latitudes).  The  model  grid  consists  of  43  equally-spaced 
longitudes  (7.  5°  spacing)  and  40  Gaussian  latitudes  ('  4 . 4°  spacing).  Thus  for  the 
longitudinal  transformation  wc  simply  use  every  third  data  point.  For  the  merid¬ 
ional  transformation  we  use  a  two-pomt  linear  interpolation  in  u  -  sin  o. 

The  second,  and  more  involved,  step  of  the  preprocessing  consists  of  inter¬ 
polating  the  data  from  the  12  mandatory  pressure  levels  to  the  model's  a  layers. 

The  analyzed  data  include  geopotential,  temperature,  wind,  and  relative  humidity 
at  all  pressure  levels  as  well  as  geopotential  at  the  surface  of  the  earth.  The  re¬ 
quired  preprocessed  data  must  consist  of  temperature,  wind,  and  specific  humid¬ 
ity  for  all  o  layers  as  well  as  pressure  at  the  surface  of  the  earth. 

We  begin  by  defining  the  model's  o -structure  (see  discussion  in  Section  2.2.2). 
Note  that  since  u  is  an  independent  variable,  can  be  defined  solely  as  a  function 
of  and  the  prescribed  u  values  at  the  bounding  surfaces  of  layer  k.  On 

the  other  hand,  p^.,  the  pressure  for  layer  k,  is  a  dependent  variable,  and  is  given 

Pj.(X,  <t>,  t)  =  p  (A,  o,  t)  Oj.  .  (72) 

In  order  to  determine  p^.(A,0,t)  we  must  first  have  p...  (A,  <J,  t).  For  the  remainder 
of  this  section,  we  shall  omit  any  reference  to  (A,q,t)  in  the  interest  of  brevity. 

W'e  determine  p..  by  applying  Eq.  (7  0)  at  the  earth's  surface,  that  is,  let 
/.  =  In  p.  and  write  Eq.  (7o)  as 

A 

|  K( Z,  -  Z )”  +  A ( Z .  -  Z)  +  ‘HZ)  -  4>,  =  0  ,  (73) 
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where  A,  B,  and  Z  are  given  by  Eq.  (71)  and 

4>(Z)  =  (*a  +  4>b)/2  -  B(Z...  -  Z)2/8  .  (74) 

Here  the  subscripts  a  and  b  indicate  values  at  the  two  mandatory  pressures  that 
are  nearest  to  the  earth's  surface.  We  can  now  easily  solve  Eq.  (73),  which  is  a 
quadratic  equation  for  Z#.  Of  the  two  solutions,  however,  only  one  is  physically 
consistent.  Once  Z^  is  known,  p^,  pk  =  p*.  <rk  and  pk  =  p*  may  then  be  computed. 
With  defined  and  p.,.  estimated,  we  may  now  proceed  to  compute  layer  values  for 
the  model  variables. 

As  mentioned  earlier,  layer  values  for  the  model  variables  are  computed  by 
interpolating  from  the  analyzed  values  at  the  mandatory  pressures.  To  do  this, 
geopotential  values  at  model  levels,  pk<  are  first  computed  by  applying  Eq.  (70) 
at  \  =  In  pk: 

$k  =  <f(Z)  +  A(Zk  -  Z)  +  I  B(Zk  -  Z)2  k  =  1.  K  ,  (75) 

where  <t>(Z)  =  (4>a  +  <J>b)/2  -  B(Zk  -  Z)2/8,  and  B,  A,  and  Z  are  given  by  Eq.  (7l). 
Here  the  subscripts  a  and  b  indicate  values  at  the  two  adjacent  mandatory  pres¬ 
sures  surrounding  pk-  Next,  layer  temperatures  Tk  are  computed  via  the  cen- 
tered-difference  form  of  Eq.  (69) 


T 


k-1 


_1 

R 


$  -  $ 
k  ^k-l 

Z.  -  Z,  . 
k  k-1 


k  =  2.  K  +  1 


(76) 


where  With  Tk  known,  layer  geopotential,  4>k,  are  obtained  by  solving 

the  finite -difference  hydrostatic  Eqs.  (46)  and  (43). 

Layer  values  for  the  two  components  of  the  pseudo-velocity,  Uk  =  uk  cos  <£ 

Vk  =  vk  cos  4>,  and  specific  humidity  Qk  are  computed  by  linear  interpolation  in  Z. 
The  appropriate  relationship  is  Eq.  (68)  where  F  represents  any  one  of  the  vari¬ 
ables  U,  V,  or  Q. 

Since  analyzed  moisture  content  is  given  in  terms  of  relative  humidity,  RH, 
we  must  convert  RH  to  Q  before  the  vertical  interpolation.  For  completeness, 
the  conversion  formulas  are  included  here: 


0.  622  e 
p  -  0.  378  e 


(77) 
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where 


e  =  e  X  RH  cb 
s 

eg  =  b  exp  (-a/T)  cb 

_  6Lv  .  2.  51  X  IQ4  o, 
a  R  1.  61  X  2. 87 


b  =  eg(273.  16)  exp  (a/273.  16)  =  0.611  exp  (a/273.  16)  cb  . 

The  preprocessing  is  completed  by  transforming  the  variables  from  model 
grid  values  to  a  set  of  truncated  spherical  harmonic  coefficients  according  to  the 
forward  transform  procedure  described  in  Section  2.2.  1.  The  only  additional 
point  to  note  here  is  that  the  spherical  harmonic  coefficients  of  vorticity  and 
divergence,  p™  and  D**1,  are  computed  spectrally  from  the  Fourier  coefficients 
of  the  pseudo-velocity  components.  This  is  done  as  follows:  the  definitions  of 
vorticity  and  divergence  art 


n  =  f  +  f 


1 

9 

a  cos'"  <t> 


av 

ax 


cos 


au ' 

3  d 


+  2  fl  sin  d 


l)  = 


1 


a  cos  d  L 


3U  ,  ,3V 

3X  +  cos  *  3d 


(78) 


and  since  these  equations  are  linear,  their  Fourier  transforms  are  simply 


in  2 

a  cos 


au 


imV  -  cos  d  -t. 


ad 


2U  sin  d  m  =  0 


m  *  0 


D 


a  cos  d 


av 


im  U  +  cos  d 


ad 


(79) 


The  spherical  harmonic  coefficients  can  be  obtained  by  applying  the  operator 


ir/2 

f  (  )  Pnl  cos  b  dd 
J  m  n 

-it/  2 
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. . .. 


to  Eq.  (79).  Noting  the  similarity  in  structure  between  Eq.  (79)  and  the  horizon¬ 
tal  advection  terms  in  the  equations  of  motion  (see  Section  2)  we  find  that  n™  and 
O™  can  be  obtained  by  use  of  the  operator  L™  defined  in  Eq.  (28),  that  is. 


m  =  Lm(v  w  j  + 

n  n  in*  m 


S-  (m,  n)  =  (0,  l) 


(m,  n)  *  (0,  1) 


Dm  =  Lm(U  ,V  ) 
n  n  in  m 


4.1.2  POSTPROCESSING 


The  problem  of  postprocessing  is  essentially  the  inverse  of  that  of  preprocess¬ 
ing.  Here  the  spherical  harmonic  coefficients  of  geopotential  and  pressure  at  the 
earth's  surface,  and  layer  values  of  temperature,  wind,  and  specific  humidity  are 
known.  The  problem  is  to  find  their  values  as  well  as  those  of  the  geopotential  at 
mandatory  levels  on  a  latitude -longitude  grid.  Briefly,  the  procedure  is  as 
follows. 

First,  the  grid  point  values  of  the  variables  on  the  o  layers  are  computed 
from  the  spectral  coefficients  by  using  the  inverse  transform  procedure  described 
in  Section  2.2.  1.  This  can  be  done  for  any  latitude -longitude  grid.  In  this  case, 
we  use  the  analysis  grid,  that  is,  a  regular  spacing  of  2.  5°  in  both  latitude  and 
longitude. 

Next  we  must  interpolate  from  the  o  layers  back  to  the  mandatory  pressure 
surfaces.  To  find  the  geopotential  at  the  mandatory  pressures,  we  follow  a  three- 
step  procedure: 

(a)  are  computed  from  T^  using  the  disc  rete  hydrostatic  equation 

VrVRTkArV  k  =  K+i,2  ,  (8i) 


where 


Zk  =  in  \  =  ln  (p*  °k* 


<p  s  <J> 

K+l  ■ 


PK+i  =  P* 


(b)  As  a  consequence  of  the  hydrostatic  approximation  and  the  assumption 
that  the  geopotential  is  piece-wise  quadratic  in  Z,  T^  is  now  given  by 


k  =  2.  K 


(82) 


T 


k 


(Tk.i  ♦  Tkl/2  ♦ 


Here  Tj,  T^+1  must  be  computed  by  extrapolation, 

(c)  With  and  known,  geopotential  values  at  mandatory  pressures 

$  ,  are  computed  by 

mand  *  J 


mand 


=  $(Z)  +  A(Z  .  -  Z)  +  i  B(Z  .  -  Z)2 
mand  2  mand 


(83) 


where  $(Z)  =  ($a  +  $b>/2  -  ®^mancj  '  Z)2/ 8  and  B,  A,  Z  are  given  as  in  Eq.  (71) 
except  that  p&,  p^  are  now  pressures  at  two  adjacent  model  interfaces  between 
which  Pman(j  lies.  In  the  case  where  Pmancj  K  P*»  we  simply  adopt  NMC’s  Shuell 
method  to  obtain  *mancj  values  at  mandatory  pressures  below  the  earth's  surface. 

For  values  of  wind,  temperature,  and  specific  humidity  at  mandatory  sur¬ 
faces,  linear  interpolation/extrapolation  in  Z  is  used: 


F 


mand 


^mand 


-  V 


(84) 


Here  Pmancj  lies  between  layer  pressures  pa  and  pfa,  and  F  represents  T,  V,  V, 

or  Q.  For  specific  humidity  we  compute  Qmancj  only  for  cases  where 

Pmand  >  300  mb"  Finally*  Qmanc|  are  to  be  converted  to  RHmanci  using  Eq.  (77). 

4.2  Initialization  Procedure 

13 

The  normal-mode  initialization  of  Machenhauer  is  applied  to  the  initial  data 

14 

following  the  procedure  developed  by  Ballish.  The  theoretical  background  and 

15 

specific  details  of  the  procedure  are  described  by  Sela. 

We  shall  confine  ourselves  to  presenting  data  that  show  the  differences  in  the 
basis  functions  between  our  model  and  the  NMC  model.  These  variations  are  due 
to  the  differences  in  sigma  structure  and  horizontal  truncation. 


13.  Machenhauer,  B.  (1977)  On  the  dynamics  of  gravity  oscillations  in  a  shallow 

water  model  with  applications  to  normal  mode  initialization,  Beit.  Phys. 
Atmos.  50:253-27 1. 

14.  Ballish,  A.  B.  (1980)  Initialization  Theory  and  Application  to  the  NMC  Spectral 

Model,  PhD  Thesis,  Dept,  of  Meteorology,  Univ.  of  Maryland. 

15.  Sela,  J.G.  (1982)  The  NMC  Spectral  Model,  NOAA  Technical  Report  NWS-30, 

3R  pp. 
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First  of  all,  as  Ballish  remarked,  there  is  no  a  priori  reason  that  the  eigen¬ 
values  associated  with  the  vertical  structure  of  any  model  be  real.  However,  their 
values  must  be  real  in  order  that  the  normal-mode  representation  be  meaningful. 
Our  experience  with  a  number  of  samples  shows  that  both  the  eigenvalues  and  the 
eigenvectors  are  robust  in  the  sense  that  they  remain  real  and  stable  with  regard 
to  changes  in  either  the  sigma  structure  or  the  basic  temperature  profile  of  the 
model.  Tables  2  and  3  present  the  sets  of  eigenvalues  for  the  six-layer  and  12- 
layer  models,  respectively,  with  three  different  basic  temperature  profiles. 
Figures  2  through  5  show  the  eigenvectors  corresponding  to  the  four  largest  eigen¬ 
values  of  some  of  the  cases  included  in  Tables  2  and  3.  Variations  of  the  eigen¬ 
values  with  structure  and  profile,  and  similarities  in  the  profiles  of  the  eigenvec¬ 
tors  among  the  various  structures  and  profiles  clearly  illustrate  these  points. 

With  the  rhomboidal  truncation  M  =  15,  the  matrix  of  the  eigensystem  for 
either  symmetric  or  antisymmetric  flow  is  of  order  24  (=3  X  8).  The  correspond¬ 
ing  eigenvectors  consist  of  eight  components  each  of  divergence,  vorticity,  and 
the  composite  function,  W  =  <t>  +  RT()q.  For  a  given  zonal  wavenumber,  both  the 
symmetric  and  antisymmetric  flows  will  generally  have  24  distinct  eigenvalues 
and  eigenvectors.  The  only  exception  is  found  with  m  =  0,  where  the  symmetric 
flow  includes  14  gravity  modes  and  two  stationary  modes,  while  the  antisymmetric 
flow  contains  16  gravity  modes.  All  Rossby  modes  are  degenerate  in  this  case. 

For  m  *  0,  of  the  24  modes,  16  are  gravity  modes  that  are  equally  divided  be¬ 
tween  eastward-  and  westward-propagating  waves  while  the  remaining  eight  are 
the  Rossby  modes.  Faster  gravity  waves  have  smaller  length  scale  in  the  merid¬ 
ional  direction,  while  faster  Rossby  waves  are  associated  with  larger  length  scale. 
Table  4  summarizes  some  of  the  characteristics  of  the  horizontal  normal  modes 
associated  with  the  M  =  15,  six-layer  AFGL  model  with  a  300°  K-isothermal  basic 
temperature  profile. 

For  the  experiments  reported  in  Section  5.2,  the  initialization  was  applied 
only  to  the  two  largest  vertical  modes  and  only  to  those  horizontal  modes  whose 
periods  are  less  than  48  hr.  The  best  results  were  obtained  with  two  iterations. 
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Table  2.  Eigenvalues  of  Vertical  Modes  of  Six -layer  Models  (m  sec  ) 


■ 

GL  Six-layer 

Equi 

-thickness  Six-layer 

1 

300° 

Isothermal 

Standard 

270° 

Isothermal 

300° 

Isothermal 

1 

93768 

104865 

116516 

93512 

104395 

115994 

2 

19945 

28085 

31206 

14589 

21193 

23547 

3 

1356 

3463 

3848 

1362 

3616 

4017 

4 

289 

743 

826 

273 

830 

923 

5 

99 

280 

311 

72 

219 

144 

6 

15 

43 

48 

12 

36 

40 
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Table  3.  Eigenvalues  of  Vertical  Modes  of  12 -layer  Models  (m  sec  ) 


1 

NMC  12 -layer 

Equi-thickness  12- 

layer 

300° 

Isothermal 

Standard 

27  0° 

Isothermal 

300° 

Isothermal 

1 

94463 

106102 

1I789I 

94140 

105484 

11 7204 

2 

27  094 

33855 

37616 

22349 

28800 

32  000 

3 

5209 

9372 

10413 

3051 

6810 

7  567 

4 

1559 

3194 

3549 

1000 

2174 

2415 

5 

696 

1470 

1634 

361 

912 

1013 

6 

324 

715 

794 

160 

426 

474 

7 

152 

388 

431 

76 

218 

242 

8 

86 

209 

233 

39 

112 

125 

9 

44 

105 

116 

19 

57 

64 

10 

22 

59 

66 

9 

27 

30 

11 

7 

21 

24 

3 

10 

11 

12 

2 

5 

6 

1 

2 

2 
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I' igurc  4.  Vertical  Profile  of  the  Eigenvectors  for 
the  NMC,  Twelve-layer,  300  K-Isothermal  Struc¬ 
ture 


Figure  5.  Vertical  Profile  of  the  Eigenvectors  for 
an  Equi-thickness,  Twelve-layer,  Standard  Struc- 


Table  4.  Some  Characteristics  of  the  Horizontal  Normal  Modes  With  16-wave 
Rhomboidal,  AFGL  Six-layer  Model,  and  the  300°-K  Isothermal  Basic  Temper¬ 
ature  Profile 


Period  in  Hours  of 

Vertical 

Zonal 

Fastest 

Slowest 

Fastest 

Mode 

Wave 

Gravity 

Gravity 

Rossby 

Number 

Number 

Mode 

Mode 

Mode 

0 

2.  09 

18.  0 

uC 

1 

1.  97 

29.  5 

27.  8 

5 

1.  58 

6.  07 

72.9 

15 

1.  07 

2.  11 

19.2 

0 

4.00 

26.  6 

jC 

1 

3.  76 

38.4 

31.  8 

5 

3.  04 

11.8 

7  5.0 

15 

2.  05 

4.  09 

193 

0 

9.  36 

46.  0 

oC 

1 

9.  00 

78.  7 

43.  0 

5 

7.81 

33.  3 

86.  3 

15 

5.  67 

11.7 

196 

0 

11.  7 

67.  9 

jC 

1 

11.  1 

125 

101 

5 

11.2 

69.  5 

104 

15 

10.  1 

25.4 

205 

0 

12.  1 

88.2 

1 

jC 

1 

11.2 

163 

137 

5 

12.2 

109 

12  0 

15 

12.5 

41.3 

216 

J 

0 

12.  3 

177 

JC 

1 

11.2 

242 

289 

5 

12.7 

246 

179 

_ _ 

14.4 

104 

2  5  7 

5.  RESULTS 

In  this  section  we  present  the  results  of  several  time  integrations  to  establish 
the  performance  of  the  low  resolution  baseline  model.  In  all  cases  we  used  a 
rhomboidal  truncation  of  M  =  15  with  six  unequally-spaced  c  layers  (see  Figure  1 
and  Table  l).  Unless  otherwise  noted,  we  used  the  semi-implicit  time  scheme 
with  a  time  step  at  4t  =  60  min.  Two  types  of  initial  data  were  used  for  testing 
purposes:  (a)  idealized  analytic  Rossby-Haurwitz  waves  and  (b)  global  l'GGE  level 
III-A  data  for  the  weeks  of  15  January  1978  and  16  July  1978. 
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.'i  l  Kiirnasls  Willi  Vnah tic  Initial  C.omlitioiM 

For  the  prelim inary  tost  of  the  model's  coding  and  performance,  we  used  the 

analytic  initial  conditions  suggested  by  Phillipslr  for  a  barotropic  model  and  mod- 

7 

ificd  for  use  in  a  multilayer  model  as  shown  by  Bourke.  The  model  was  initialized 
with  the  same  flow  pattern  in  all  layers.  It  consisted  of  a  mean  zonal  flow  plus  a 
zonal  wavenumber  4  eddy  with  the  streamfunetion  given  by 

e(X,  o)  -  C?  P?(sin  o)  +  <i5  pt  (sin  d)  e4^  +  i/l  p1  (sin  <?)  ,  (85) 

II  J  )  J  o 

where  the  spectral  amplitudes  are 

0  2 
c  j  =  -a  u' 


4  .  a2u  /  2  •  r>! 

fc.)  ‘  1890  V  ii 


t ) 

with  .  -  ~  nncj  is  the  angular  velocity  of  the  earth.  Additional  specifications 
of  the  initial  conditions  included  an  isothermal  atmosphere  with  T  =  266.  5°K  in 
all  layers,  a  mean  surface  pressure  of  100  cb,  and  no  topography  so  that  -  0 
everywhere.  The  initialization  was  completed  with  a  non-divergent  balance  con¬ 
dition  obtained  from  tlq.  (24)  by  initially  setting 


_ n_ 

at 


4>m  =  i) 


and  then  solving  for  q™.  In  this  ease,  the  model  was  treated  as  inviseid,  adiabatic, 
and  dry  (that  is,  all  of  the  physical  parameterizations  were  turned  off). 

We  note  that  the  analytic  barotropic  non-divergent  solution  consists  of  a  zonal 
wavenumber  4  pattern  that  preserves  its  shape  and  amplitude  and  rotates  from 
west  to  east  with  a  phase  speed  of  9.6°  per  day.  For  the  multilayer  model,  the 
flow  should  resemble  the  barotropic  solution  for  the  first  few  time-steps  and  then 
slowly  deviate  as  divergence  and  vertical  shear  develop. 

4 

In  Table  5,  we  show  the  amplitude  and  phase  speed  of  the  ^  mode  resulting 
from  two  different  72 -hr  forecasts,  one  using  an  explicit  time  scheme  with  a  time- 
step  of  10  min  and  the  other  using  the  semi -implicit  scheme  described  before. 


16.  Phillips,  N.A.  (1959)  Numerical  integrations  of  the  primitive  equations  on  the 
hemisphere,  Mon.  Wea.  Rev.  87:333-345. 
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Table  5.  Amplitude  (10*  in"  sec  *)  and  Phase  Speed  (Degrees  per  Day)  of 
the  if..  Mode 


Kxplicit  Time  Scheme  (At 

=  10  min) 

24  hr 

48  hr 

72 

hr 

a 

A 

P 

A 

P 

A 

P 

0.  :>022 

3.  995 

9.  47 

3.  949 

9.  48 

3.  838 

9.  49 

II.  IPS') 

3.  999 

9.  37 

3.  87  9 

9.  30 

3.  826 

9.  15 

0.  36PP 

3.  986 

9.  13 

3.  951 

8.  97 

3.951 

8.  f7 

0.  622  0 

4.  016 

8.  94 

4.011 

8.  90 

4.  044 

8.  92 

o.  8242 

4.  030 

8.  85 

4.  04  a 

8.  91 

4.  041 

8.  95 

0.  94  97 

4.  639 

8.  84 

4.  029 

8.  95 

3.  983 

8.  93 

Average 

4.  009 

9.  09 

3.  981 

9.  05 

3.  961 

9.  05 

Analytic 

4.  023 

9.  6 

4.  023 

9.  6 

4.  023 

9.  6 

Semi-implicit  Time  Scheme  (At  =  60  min) 

24  hr 

48  hr 

72  hr 

a 

A 

P 

A 

P 

A  1 

P 

0.  0622 

4.  001 

9.  46 

3.  957 

9.41 

3.825 

9.  37 

(I.  1985 

4.  003 

9.  37 

3.  901 

9.  20 

3.  854 

9.  00 

0.  3699 

3.  992 

9.  10 

3.  975 

8.  86 

3.  996 

8.78 

0.  6220 

4.  030 

8.  93 

4.  040 

8.  83 

4.  113 

8.75 

0.  8242 

4.  051 

8.  85 

4.091 

8.  85 

4.  079 

8.  82 

0.  9497 

4.  063 

8.  85 

4.  093 

8.  88 

4.  038 

8.  80 

Average 

4.  020 

9.  08 

4.  010 

8.  97 

4.  002 

8.89 

In  both  cases,  we  find  that  the  model  behaves  as  expected,  that  is,  the  flow  re¬ 
sembles  the  barotropic  solution  with  deviations  that  increase  with  time.  As  noted 
7 

by  Boorke,  we  also  find  that  the  upper  layers  tend  to  preserve  the  phase  speed  of 

the  barotropic  flow  while  the  lower  layers  tend  to  preserve  the  amplitude  of  the 

wave.  In  general,  the  vertically-averaged  amplitude  and  phase  speed  both  decrease 

slightly  as  the  forecasts  proceed.  The  small  differences  between  the  explicity  and 

semi-implicit  results  can  easily  be  explained  in  terms  of  the  known  characteris- 

17 

tics  of  each  of  these  time  schemes  (for  example,  Kurihara  ). 


17.  Kurihara,  Y.  (1965)  On  the  use  of  implicit  and  iterative  methods  for  the  time 
integration  of  the  wave  equation,  Mon.  Wea.  Rev.  93:33-46. 
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Finally,  the  analytic  initialization  experiments  provide  an  opportunity  to 
check  the  conservation  properties  of  the  model.  For  the  72 -hr  semi-implicit 
integration,  the  changes  in  total  energy  (sum  of  kinetic,  internal,  and  gravitational 
potential),  mass,  and  angular  momentum  were  -0.0034,  -0.0035,  and  -0.  162  per¬ 
cent,  respectively.  For  the  explicit  integration  the  change  in  energy  was  -0.  0026 
percent.  All  of  these  values  are  considered  acceptable  and  allow  us  to  proceed  to 
the  more  realistic  simulations  that  include  all  of  the  physical  parameterizations 
and  use  actual  data  for  initial  conditions. 


.1.2  Forecasts  With  Global  Data 

As  was  mentioned  previously,  the  data  used  for  the  global  test  forecasts  con¬ 
sisted  of  FGGE  level  1II-A  data  sets  for  the  weeks  beginning  15  January  1078  and 
16  July  1978.  In  view  of  the  limited  number  of  days  of  data  available,  the  decision 
was  made  to  run  four  test  forecasts  —  one  for  the  first  48-hr  period  of  each  week 
and  one  for  the  last  48-hr  period  of  each  week.  This  was  done  in  order  to  test  the 
model  on  as  many  synopticaliy  independent  cases  as  possible.  While  the  separa¬ 
tion  of  four  days  may  not  be  quite  enough  for  synoptic  independence,  it  was  the  best 
that  could  be  done  under  the  circumstances.  Thus  the  four  test  forecasts  included 
two  winter  cases  with  00Z  15  January  1978  and  12Z  19  January  1978  and  two 
summer  cases  with  DoZ  16  July  1978  and  12Z  20  July  1978  as  the  initial  times. 

We  note  that  in  the  six -layer,  M  =  15  rhomboidal  configuration  with  At  =  60  min 
the  model  required  83K  words  of  memory  and  used  18  min  of  GPU  time  for  each 
24  hr  of  simulation  on  the  COG  6600. 

In  the  course  of  evaluating  the  model's  performance,  we  found  it  helpful  and 
interesting  to  assess  the  relative  influences  of  each  of  the  three  components  of  the 
model  (preprocessing,  prognostication,  and  postprocessing)  as  shown  in  Figure  6. 
It  is  obvious  that  the  preprocessing  and  postprocessing  are  purely  static  while  the 
prognostication  includes  both  static  and  dynamic  changes.  If  the  errors  accrued 
in  the  static  transformation  during  the  processing  stages  were  statistically  inde¬ 
pendent  of  those  produced  during  the  prognostication  the  total  error  variance  of 
the  model  would  be  the  sum  of  the  two  partial  error  variances.  In  practice,  how¬ 
ever,  because  of  the  presence  of  common  elements  in  both  physical  assumptions 
and  mathematical  procedures  among  different  components  we  do  not  expect  a 
strict  statistical  independence  between  the  processing  parts  and  the  prognostic 
part.  Nonetheless,  it  was  thought  that  such  a  partition  might  shed  a  light  into  the 
composition  and  nature  of  the  model  performance.  Toward  this  end,  the  perform¬ 
ance  was  measured  by  the  RMSD  between  the  model-produced  values  and  the 
corresponding  analyzed  or  "truth"  values  at  the  grid  points.  Thus, 
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(A) 


(S) 


Figure  fi.  Schematic  Diagram  of  Model 
Evaluation 


KMS1) 


where  My  and  T„  are  the  model-produced  and  "truth"  values  at  grid  point  (i,  1)  in 
which  i  =  l,  2,  . . .  ,  1  and  j  s  1,  2,  ,,,,  J  are  the  longitudinal  and  latitudinal 
indices,  respectively,  and  V\\  is  the  Gaussian  weight  at  latitude  i. 

For  the  total  performance  of  the  model,  the  values  of  M  were  taken  from 
the  model  forecast  and  TV  were  the  values  interpolated  from  the  FGGE  III -A  data 
set  to  the  Gaussian  grid,  hereafter  referred  to  as  the  analyzed  values.  For  the 
performance  of  the  processing  parts.  My  were  synthesized  from  the  analyzed 
values  at  the  time  of  verification  by  putting  the  latter  through  the  preprocessing 
and  postprocessing  parts  while  bypassing  the  prognostication  of  the  model,  i  inallv, 
the  performance  of  the  prognostication  is  represented  by  the  KMSD  between  the 
model-forecast  values  and  the  synthesized  values.  A  table  of  global  performances 
for  vector  wind,  temperature,  and  height  at  three  different  pressure  levels  sum¬ 
marizes  the  model  evaluation  in  each  of  the  four  test  runs  (Table  c  through  Table  "). 
In  these  tables  the  letters  A,  S,  F  refer  to  the  analyzed,  synthesized,  and  fore¬ 
casted  fields,  respectively. 


■  ,  ■  *  4*  -  ,  .  y. 
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Table  6.  Performance  (RMS  Differences) 


15  January  1978,  00Z 

Level 

V(m  sec  *) 

T  (°K> 

Z  (m) 

(mb) 

Pair 

24  hr 

48  hr 

24  hr 

48  hr 

24  hr 

48  hr 

9 

5,4 

6.  8 

2.  6 

3.  6 

31.  6 

42.2 

850 

2.  5 

2.  5 

1.7 

10.4 

4.  8 

6.4 

2.4 

41.  7 

6.7 

9.  0 

1.  9 

57.  3 

500 

3.  6 

3.7 

n 

12. 7 

Plllll 

5.  5 

7.9 

m 

56.  9 

EBi 

11.3 

4.  4 

4.  8 

58.  8 

88.7 

250 

I 

5.  5 

2.8 

2.6 

14.2 

13.  6 

8.6 

12.2 

2.3 

_ 

3.  0 

56.  5 

_ 

86.  8 

Table  7.  Performance  (RMS  Differences) 


19  January  1978,  12  Z 


Level 

(mb) 

Pair 

V  (m  sec  *) 

T  (°K) 

Z  (m) 

24  hr 

48  hr 

24  hr 

48  hr 

24  hr 

48  hr 

(F,A) 

5.3 

6.7 

2.5 

3.4 

850 

(S,  A) 

2.  5 

2.8 

1.6 

1.  9 

(F,S) 

4.  6 

6.  0 

2.2 

3.  3 

7.  1 

8.7 

2.0 

2.6 

38.4 

50.7 

500 

H 

3.7 

3.7 

in 

1.  5 

10.  4 

12.  0 

EHI 

5.8 

7.  6 

2.3 

38.  0 

50.2 

B 

11.7 

13.  9 

4.3 

4.  8 

57.  1 

81.  8 

250 

K 

5.  9 

5.  5 

2.  8 

2.8 

14. 0 

13.  5 

[pill 

8.  8 

11.6 

2.2 

2.8 

54.8 

79.  8 
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Table  8.  Performance  (RMS  Differences) 


16  July  1978,  002 


JLevel 

V  (m 

sec  *) 

T 

— 

°K) 

2 

(m) 

(mb) 

Pair 

24  hr 

48  hr 

24  hr 

48  hr 

24 

hr 

48  hr 

| 

4.  7 

6.  8 

2.3 

25. 

7 

41.7 

850 

V 

2.3 

1.6 

9. 

6 

8.6 

4.  1 

2.0 

2.  9 

25. 

3 

41.2 

(F.A) 

8.  6 

1 

34. 

7 

55.4 

500 

(S,  A) 

3.4 

11. 

1 

11.  0 

<F,S) 

7.  8 

34. 

0 

55.  0 

(F.A) 

14.2 

4.4 

4.7 

55. 

4 

83.8 

250 

(S.A) 

6.  3 

3.0 

3.  0 

15. 

0 

15.  5 

<F,S) 

2.2 

2.  6 

53. 

0 

81.  9 

Table  9.  Performance  (RMS  Differences) 


20  July  1978,  122 

Level 

_7 

V  (m  sec  ) 

T  (°K) 

2(m> 

(mb) 

Pair 

24  hr 

48  hr 

24  hr 

48  hr 

24  hr 

48  hr 

(F.A) 

5.0 

6.  3 

2.7 

3.  1 

■ 

37.7 

850 

(S.A) 

2.  3 

2.  1 

1.5 

1.  5 

8.  8 

(F.S) 

4.  5 

5.9 

2.  5 

3.  1 

37.  1 

(F.A) 

6.2 

8.  1 

2.  1 

52.2 

500 

(S.A) 

3.  3 

3.4 

1.6 

n 

9.8 

(F.S) 

5.2 

7.  1 

1.  6 

52.  1 

(F.A) 

11.  0 

12. 8 

4.4 

4.8 

59.4 

86.  1 

250 

(S.A) 

6.2 

6.2 

3.  0 

3.  0 

15.  3 

15.2 

(F.S) 

7.9 

10.4 

2.2 

_ 

2.5 

_ 

57.4 

84.  5 

Three  features  common  to  all  runs  seem  to  stand  out  well.  First,  the  total 
MSD  are  approximately  equal  to  the  sum  of  that  due  to  the  processings  and  that 
due  to  the  prognostication.  Second,  the  performance  of  the  processing  is  quali¬ 
tatively  invariant  with  respect  to  time  and  case,  indicating  that  the  measure  at  the 
global  level  is  insensitive  to  the  synoptic  details.  The  last  feature  that  follows  as 
a  combined  result  of  the  previous  two  characteristics  is  that  the  degradation  with 
time  of  the  total  model  performance  is  due  mostly  to  the  prognostication.  The 
values  of  performance  are  comparable  in  magnitudes  with  the  statistics  of  other 
current  operational  and  research"  models. 

For  a  qualitative  assessment  of  the  model's  performance,  in  Figure  7  we 
show  the  analyzed  1000-  and  5U0-mb  heights  for  00Z  15  January-17  January  and 
the  corresponding  24-  and  48-hr  forecasts  (16  January  and  17  January,  respec¬ 
tively)  for  the  Northern  Hemisphere.  In  general,  the  results  are  as  might  be 

expected  based  on  model  performance  reported  by  other  researchers  using  com* 
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parable  models.  ’  *  At  24  hr  (Figures  7b  and  7c)  the  model  is  able  to  correctly 
forecast  the  location  and  movement  of  the  major  surface  systems  and  the  500  mb 
troughs  and  ridges.  At  48  hr  (Figures  7d  and  7e)  the  global  patterns  are  predicted 
reasonably  well  but  phase  errors  start  to  become  noticeable.  Also,  a  major 
problem  of  a  low  resolution  model  is  its  inability  to  accurately  treat  smaller  scale 
features  such  as  cutoff  lows  and  sharply  defined  fronts. 

6.  CONTUSION 

The  formulation  of  the  baseline  version  of  the  AFGL  moist  global  model  has 
been  described  in  detail.  It  is  patterned  after  the  NMC  spectral  model.  *  The 
numerical  methods  in  the  model  include  spectral  representation  in  the  horizontal, 
the  Arakawa  vertical  differencing  scheme,  and  the  semi-implicit  time  scheme. 

The  data  processing,  normal -mode  initialization,  and  the  parameterized  physical 
processes  have  all  been  adapted  from  the  NMC  model.  Test  forecasts  for  a  limited 
number  of  cases  indicate  that  the  model  performs  as  well  as  other  comparable 
large-scale  models. 

6.1  Future  Research 

It  should  be  emphasized  that  the  model  documented  in  this  report  is  a  baseline 
version  and  is  therefore  only  a  starting  point  for  additional  research.  One  of  the 
major  goals  of  our  future  effort  is  to  provide  the  model  with  a  useful  cloud 

NMC  global  spectral  model  and  GWC  AWSPE,  private  communication. 
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OBSERVED  HEIGHT  FIELD 
OOZ  500mb 

15  JANUARY  1978 
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OBSERVED  HEIGHT  FIELD 
OOZ  lOOOmb 

15  JANUARY  1978 

Figure  7a.  Analyzed  Height  Fields  for 
1000  mb  and  500  mb  for  OOZ  15  January 
1978  (Initial  Fields) 
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OBSERVED  HEIGHT  FIELD 
00 Z  1 000 mb 

16  JANUARY  1978 


Figure  7b.  Analyzed  Height  Fields  for 
1000  mb  and  500  mb  for  00Z  ]6  January 


OOZ  500 mb 

16  JANUARY  1978 


FORECAST  HEIGHT  FIELD 
OOZ  lOOOmb 

16  JANUARY  1978 


Figure  7c.  Height  Fields  for  24 -hr  Forecast 
for  00 Z  16  January  1978 
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FORECAST  HEIGHT  FIELD 
OOZ  500  mb 

17  JANUARY  1978 
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FORECAST  HEIGHT  FIELD 
OOZ  lOOOmb 

17  JANUARY  1978 


Figure  7e.  Height  Fields  for  48-hr  Forecast 
for  OOZ  17  January  1978 


forecasting  capability.  Some  of  the  preliminary  steps  and  procedures  necessary 

1 8 

for  this  have  already  been  described  by  Mitchell. 

Success  in  this  area  will  require  major  changes  and  improvements  in  all 
aspects  of  the  baseline  model: 

(a)  In  the  area  of  dynamics  anti  numerical  methods,  we  are  currently  working 
on  expanding  the  resolution  of  the  model  (at  least  12  layers  with  M  =  30)  as  well  as 
vectorizing  the  code  for  efficient  processing  on  the  CRAY-1.  We  are  also  consider 
ing  alternative  and  potentially  more  accurate  numerical  solutions  of  the  hydro¬ 
static  equation. 

(b)  Data  processing,  including  analysis  and  initialization,  is  one  area  that  is 
currently  receiving  substantial  attention  throughout  the  meteorological  community. 
For  our  model,  an  effort  is  under  way  to  develop  a  scheme  that  combines  optimal 
interpolation  analysis  and  nonlinear  normal  mode  initialization.  Other  research 
in  the  processing  procedures  will  focus  on  improving  the  pressure-to-sigma  and 
sigma-to-pressure  interpolation  schemes  to  make  them  more  accurate  and  con¬ 
sistent  with  the  vertical  structure  of  the  forecast  model. 

(c)  The  parameterized  physics  in  the  model  will  also  be  completely  changed. 
Efforts  are  currently  in  progress  to  develop  improved  schemes  for  boundary- 
layer  and  moisture  physics.  In  addition,  a  radiation  package  will  be  added  to  the 
model. 

(d)  Finally,  an  effort  is  being  made  to  develop  preliminary  software  that  will 
produce  derived  cloud  forecasts  based  on  the  model-forecasted  moisture  field. 


18.  Mitchell,  K.E.  (1982)  Cloud  Forecast  Fields  Comparison  Test.  AFGWC 
Technical  Note  82/ooS,  pp. 
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Appendix  A 

Numerical  Values  of  Physical  Constants 


g 

a 

il 

R  for  dry  air 

C  for  dry  air 
p 

e  ratio  of  molecular 
weight  of  water  to 
that  of  air 

R  for  water  vapor 

L  of  condensation 

L  ol  sublimation 

specific  heat  of  liquid 
water 

specific  heat  of  ice 

C  for  water  vapor 
P 


!>.  8  m  s 
fi.  371  x  iofl  m 
7.292116  X  10'3  s'1 
287. 05  J  kg’1  K"1 
1005  J  kg"1  K'1 

0.  622 

461.  5  J  kg'1  K'1 
2.  51  X  106  J  kg"1 
2.835  X  10B  J  kg"1 

4187  J  kg"1  K'1 
1922  J  kg"1  K'1 
1876.  5  J  kg"1  K"1 
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j 

J 


i 


A 
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GO 


T 


We  now  let 


T  =  aT  (1  -  a  )T 


D  =  aD  (l  -  a)D‘ 


q  =  aq  (1  -  a)q 


(B6) 


for  the  undifferentiated  linear  terms  in  Eqs.  <B2),  (B3),  and  (B5)  and  we  approx¬ 
imate  the  time  derivatives  by  centered-differencing. 


3F  F  -  F 


3t  2  At 


(B7 ) 


Equations  (B5),  (B2),  and  (B3)  may  then  be  written  as 


D+  -  2  At  a  <AT+  +  Rq+TQ)  =  D*  +  2  At  jc  +  -(n2+-^  e 


+  2At(l  -  a)  P(~  ^  (AT'  +  Rq‘Tn) 
a 


(B8) 


T+  =  T*  -  2 AtaCD+  -  2At(l  -  a)CD_  +  2Aty  , 
q+  =  q  -  2AtlarD+  +  (1  -  a£rD  )  +  2Atz 


(B9) 

(BIO) 


Let  I  be  the  identity  matrix  of  order  K.  Elimination  of  T+  and  q+  in  Eq.  (B8) 
leads  to  the  linear  system 


I  +  (2Ata)2  n(n/  —  (AC  +  RT(lr) 
a2 


D+  = 


I  -  (2Ato)2  n(n-2+  1}  (AC  +  RT()  r) 


D  + 


2Atjc  +  +  A(T‘  +  2Atc^)  +  R(q‘  +  2Atoz)X0  ,  (Bll) 

& 
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List  of 


a  radius  of  the  earth 

f  Coriolis  parameter 

g  gravity 

i  V-l 

j,k,  1  indices 

k  vertical  unit  vector 

m  order  of  spherical  harmonic  (superscript) 

n  degree  of  spherical  harmonic  (subscript) 

p  pressure 

p*  surface  pressure 

q  In  P* 

t  time 

u,v  eastward  and  northward  components  of  velocity 

A,  B  nonlinear  advection  terms 

Cp  specific  heat  at  constant  pressure 

D  horizontal  divergence 

E  kinetic  energy 
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H 

n 

I 

J 

K 

L 

v 

M 

N 


0 

u.v 
w. 

J 

Yr 

r 

? 
n 
0 


J 

fm 

n 


A 


o 

6 

e 


X 

<l> 


diabatic  heating 

meridional  derivative  of  the  associated  Legendre  function 

nonlinear  term  in  the  thermodynamic  equation 

nonlinear  term  in  the  moisture  equation 

number  of  vertical  layers 

latent  heat  of  condensation  of  water  vapor 

truncation  value  of  m 

truncation  value  of  n 

associated  Legendre  function  of  order  m  and  degree  n 

specific  humidity 

gas  constant  for  dry  air 

temperature 

basic  state  temperature 

pseudo-velocity  =  u  cos  <f>,  v  cos  <t> 

Gaussian  weight 

spherical  harmonic  of  order  m  and  degree  n 

relative  vorticity 

absolute  vorticity 

potential  temperature 

R/C 

P 

longitude 
sin  <t> 

diffusion  coefficient 
vertical  coordinate 
vertical  velocity  = 

ratio  of  molecular  weight  of  water  vapor  to  that  of  air 
latitude 

velocity  potential 
streamfunction 


64 


<f> 

<  c 

(T 

(  >k 
L) 


geopotential 

surface  geopotential 

angular  velocity  of  the  earth's  rotation 

spherical  harmonic  coefficient  of  order  m  and  degree  n 

vertical  integral 

layer  value 

level  (interface)  value 


vector 
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